library(stars) ; library(sf)
library(dplyr) ; library(tidyr)
library(ggplot2) ; library(gganimate)
library(lubridate) ; library(stringr)
This document demonstrates how to visualize spatialtemporal arrays (time-series of rasters) in R.
The example data is the daily OMI ColumnAmountNO2TropCloudScreened data, which has one .tif file for each day. There are 365 files in total.
files_OMI <- list.files("1_data/raw/OMI-NO2" , # directory of the files
pattern = "_AOI.tif$" , # filtering files with naming patterns (optional)
full.names = TRUE) # with directory path
Here’s an example of the first 6 files of the total 365 files.
head(files_OMI)
## [1] "1_data/raw/OMI-NO2/OMI-Aura_L3-OMI_MINDS_NO2d_20190101_AOI.tif"
## [2] "1_data/raw/OMI-NO2/OMI-Aura_L3-OMI_MINDS_NO2d_20190102_AOI.tif"
## [3] "1_data/raw/OMI-NO2/OMI-Aura_L3-OMI_MINDS_NO2d_20190103_AOI.tif"
## [4] "1_data/raw/OMI-NO2/OMI-Aura_L3-OMI_MINDS_NO2d_20190104_AOI.tif"
## [5] "1_data/raw/OMI-NO2/OMI-Aura_L3-OMI_MINDS_NO2d_20190105_AOI.tif"
## [6] "1_data/raw/OMI-NO2/OMI-Aura_L3-OMI_MINDS_NO2d_20190106_AOI.tif"
Since all the 365 raster files have the same spatial dimension, we can import them once and as a spatialtemporal array using the package stars.
OMI <- read_stars(files_OMI) # bad example
dim(OMI)
## x y
## 24 13
This results in a stars object with two dimensions (x,y) and 365 attributes (bands). What we want is one attribute and 3 dimensions (x, y, date), so we need to assign the along argument in read_stars as the date.
# I first extract the date from the file names (OMI-Aura_L3-OMI_MINDS_NO2d_YYYYMMDD_AOI.tif)
dates_OMI <- basename(files_OMI) %>%
str_extract("\\d{8}") %>%
as_date()
# alternatively, you can assign the time stamps directly:
# dates_OMI <- seq(as_date("2019-01-01") , as_date("2019-12-31") , "1 day")
# then assign the date values to read_stars to set the third dimension
OMI <- read_stars(files_OMI ,
along = list(date = dates_OMI)) %>% # used a named list to assign the dimension
# rename the attibute name
setNames("value")
OMI
## stars object with 3 dimensions and 1 attribute
## attribute(s):
## value
## Min. :-9.982e+14
## 1st Qu.: 4.562e+14
## Median : 1.389e+15
## Mean : 1.968e+15
## 3rd Qu.: 2.728e+15
## Max. : 5.517e+16
## NA's :83574
## dimension(s):
## from to offset delta refsys point values x/y
## x 1 24 5.25 0.25 GEOGCS["unnamed ellipse",... FALSE NULL [x]
## y 1 13 45.25 0.25 GEOGCS["unnamed ellipse",... FALSE NULL [y]
## date 1 365 2019-01-01 1 days Date NA NULL
Here OMI is a stars with three dimensions: x, y, and date.
I also import a shapefile of Switzerland to demonstrate how to add polygons on top of the rasters.
CH <- st_read("1_data/raw/Switzerland_shapefile/CHE_adm0.shp") %>%
# exclude the attribute table
st_geometry() %>%
# simplify to facilitate visualization
st_simplify(preserveTopology = TRUE , dTolerance = 0.05)
## Reading layer `CHE_adm0' from data source `/Masterarbeit/analysis/1_data/raw/Switzerland_shapefile/CHE_adm0.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 70 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 5.956063 ymin: 45.81706 xmax: 10.49511 ymax: 47.82369
## CRS: 4326
## Warning in st_simplify.sfc(., preserveTopology = TRUE, dTolerance = 0.05):
## st_simplify does not correctly simplify longitude/latitude data, dTolerance
## needs to be in decimal degrees
ggplotggplot() + geom_stars() to map stars objects in ggplot
facet_grid() to make separate plots according to dimensions (date)+ geom_sf() to add polygons# for demonstration, I only visualize a 5-day time series
selected_dates <- interval("2019-07-15" , "2019-07-19")
# I subset the original array using the date dimension value by:
# OMI %>%
# filter(date %within% selected_dates)
ggplot() +
# plot the rasters
geom_stars(
data = OMI %>%
filter(date %within% selected_dates)
) +
# plot the polygon
geom_sf(data = CH , fill = NA , color = "white") +
# one plot per date (based on the dimension date of OMI)
facet_grid(~date) +
# set color for the rasters (optional)
scale_fill_viridis_c() + # high contrast color option
# optional settings
coord_sf(expand = FALSE) +
# labels (optional)
labs(x = "Longtitude" , y = "Latitude" , fill = "OMI-NO2")
ggplot and gganimateggplot, just replace the facet_grid(~date) with transition_time(date), which tells gganimate on which dimension the animation is based.+labs(title = "Date: {frame_time})" to add the date of each frame as titlelimits argument in scale_fill...() as for example scale_fill_viridis_c(limits = c(0,1e16))# assign the plot to an object ()
OMI_animation <- ggplot() +
# plot the rasters
geom_stars(
data = OMI
) +
# plot the polygon
geom_sf(data = CH , fill = NA , color = "white") +
# one frame per date (based on the dimension date of OMI)
transition_time(date) +
# set color for the rasters (optional)
scale_fill_viridis_c() + # high contrast color option
# optional settings
coord_sf(expand = FALSE) +
# labels (optional)
labs(x = "Longtitude" , y = "Latitude" , fill = "OMI-NO2" ,
title = "Date: {frame_time}") # use this to add the date as title
Because we have 365 days, the default pace of the animation is too fast. Therefore we assign the ggplot+... to an object and manually control the pace.
animate(OMI_animation ,
duration = 365 , # 365 seconds in total
fps = 1) # one frame per second